clc;
clear;

format shortE

%% example 1 何晓明
function_u = @(x) x.*cos(x);
function_cor = @(x) exp(x);
function_f = @(x) -exp(x).*(cos(x) - 2*sin(x) - x.*cos(x) - x.*sin(x));
function_g = function_u;
function_r = function_symbolic_computing_1d(function_u, "dx");

%% example 2
% function_u = @(x) x.*x;
% function_cor = @(x) exp(x);
% function_f = function_symbolic_computing_1d(@(x) function_cor(x).*function_symbolic_computing_1d(function_u, "dx"), "dx");
% function_f = @(x) -function_f(x);
% function_g = function_u;
% function_r = function_symbolic_computing_1d(function_u, "dx");

%% example 3 华冬英
% function_u = @(x) exp(2*x).*sin(4*pi*x);
% function_cor = @(x) 1;
% function_f = @(x) exp(2*x).*((16*pi^2-4)*sin(4*pi*x) - 16*pi*cos(4*pi*x));
% function_g = @(x) 0;
% function_r = function_symbolic_computing_1d(function_u, "dx");

%% example 4 陆金甫
% function_cor = @(x) 1;
% function_f = @(x) 2;
% function_g = @(x) 0;
% function_r = @(x) 0;

%% example 5 test
% function_u = @(x) x.*x;
% function_cor = @(x) 1;
% function_f = @(x) -2;
% function_g = @(x) x.*x;
% function_r = @(x) 2.*x;

startTime = cputime;

left = 0.0;
right = 1.0;
vector_N = [8 16 32 64];

basis_type = "P1";
N_gauss_int1d = 9;
boundary_type = "Dirichlet";
% boundary_type = "NeumanLeft";
% boundary_type = "NeumanRight";
% boundary_type = "Robin";

% plot_basis_function_1d(basis_type);

error_max = zeros(size(vector_N,2), 1);
error_Linf = zeros(size(vector_N,2), 1);
error_L2 = zeros(size(vector_N,2), 1);
error_H1 = zeros(size(vector_N,2), 1);

for i = 1:size(vector_N,2)
    N = vector_N(i);

    [uh, U] = poisson_1d(left, right, N, basis_type, N_gauss_int1d, function_cor, function_f, function_g, function_r, boundary_type);

    error_max(i,1) = processing_error_1d(left, right, N, basis_type, N_gauss_int1d, function_u, uh, "Max");
    error_Linf(i,1) = processing_error_1d(left, right, N, basis_type, N_gauss_int1d, function_u, U, "Linf");
    error_L2(i,1) = processing_error_1d(left, right, N, basis_type, N_gauss_int1d, function_u, U, "L2");
    error_H1(i,1) = processing_error_1d(left, right, N, basis_type, N_gauss_int1d, function_r, U, "H1");

    fprintf("完成%d次\n", i);
end

% plot_function_1d(uh, left, right);

processing_convergence_rate_1d(error_max, vector_N);
processing_convergence_rate_1d(error_Linf, vector_N);
processing_convergence_rate_1d(error_L2, vector_N);
processing_convergence_rate_1d(error_H1, vector_N);

% plot_function_1d(function_u, node);

elapsedCPUTime = cputime - startTime;
disp(['Elapsed CPU time is ', num2str(elapsedCPUTime), ' seconds.']);